
maindir<-getwd()

modeltype<-"simplified"

source("./auxfiles/loaddata.R")
lowermargin<-1.1
uppermargin<-1.1
goldtol <- params$goldsearch_tolerance
shocks.wage<-shocks.wagereal
shocks.spwage<-shocks.spwagereal
setwd(maindir)

#sourceCpp("~/Dropbox/DynamicDI/StructuralModel/cpp/Valsim_corrwage_wpolicy.cpp", rebuild=TRUE)


temp<-as.data.table(read.csv(paste0("./output/empiricalhealthprobs_","oldsample","_marginal.csv"),row.names=NULL))
healthprob_old<-list()
  for(typ in 1:numtypes){
    healthprob_old[[typ]]<-array(NA,dim=c(50*agelength,6,6))
    for(s in 1:3){
      for(sp in 1:2){
        pos<-1
        for(t in 1:50){
          for(rep in 1:agelength){
            healthprob_old[[typ]][pos,(s-1)*2+sp,]<-as.vector(t(outer(temp[type==typ&from==s & age == t+22,p]
                                                                  ,temp[type==typ&from==sp & age == t+22 & to !=3,sp_p]))) #probabilities  sum to 1
            pos<-pos+1
          }
        }
      }
    }
  }



skip<-1
if(skip==0){
  simdata<-Valsim_solve(as.vector(param_input),params,healthprob_old, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                        shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                        marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=1)
  
  simdata<-cleansim(simdata,select=select,selectC=selectC)
  
  ###WILLINGNESS TO PAY FOR POLICY REFORMS, OLD HEALTH PROBS
  wtpold<-evalwtp(paramvals=newguess,
               marriageprob=marriageprob, healthprob=healthprob_old,
               foodstamps=foodstamps, numkids=numkids, oecd=oecd,
               Wage=Wage, spWage=spWage, select=select,
               selectC=selectC,
               marriageprob_init=marriageprob_init,
               shocks.wage = shocks.wage, shocks.spwage = shocks.spwage,
               wagepars=wagepars, spwagepars=spwagepars,
               onlyFT=1,evalasset=FALSE, onlySS=FALSE,onlymeanstest=FALSE,
               single_sameallowanceprob=single_sameallowanceprob,
               health_simpleprefs=health_simpleprefs, single_sameprefs=single_sameprefs)
  saveRDS(wtpold,paste0("./modeloutput/wtp_old_",modeltype,".RDS"))

  #WTP TO BE ON DI
}



# 
# rm(wtp)
# gc()
# 
# saveRDS(output,"./figures/wtp_decompose/wtp_decompose_oldhealthprob.RDS")
# rm(output)

#setThreadOptions(numThreads=20,stackSize=1e07)
setThreadOptions(numThreads=20,stackSize=1e07)


setwd("./modeloutput/")
simdata<-Valsim_solve(as.vector(param_input),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                      shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                      marriageprob_init,numsims,lumpsum=0,print=1,workDI=0,onlyFT=onlyFT)
setwd("..")

simdata<-cleansim(simdata,select=select,selectC=selectC)

skip<-1
if(skip==0){
  ###WILLINGNESS TO PAY FOR POLICY REFORMS
  wtp<-evalwtp(paramvals=newguess,
               marriageprob=marriageprob, healthprob=healthprob,
               foodstamps=foodstamps, numkids=numkids, oecd=oecd,
               Wage=Wage, spWage=spWage, select=select,
               selectC=selectC,
               marriageprob_init=marriageprob_init,
               shocks.wage = shocks.wage, shocks.spwage = shocks.spwage,
               wagepars=wagepars, spwagepars=spwagepars,
               onlyFT=1,evalasset=FALSE, onlySS=FALSE,onlymeanstest=FALSE,
               single_sameallowanceprob=single_sameallowanceprob,
               health_simpleprefs=health_simpleprefs, single_sameprefs=single_sameprefs,
               filename="simplified")
  
  saveRDS(wtp,paste0("./modeloutput/wtp_",modeltype,".RDS"))
  
  #WTP TO BE ON DI
} else { wtp<-readRDS(paste0("./modeloutput/wtp_",modeltype,".RDS")) }
temp <-read.csv(paste0("./V_orig/","V.csv"),header=FALSE)[,1:numtypes]
V0<-array(NA,c(numtypes,50,dim(Wage)[3],dim(spWage)[3],dim(Asset)[3],2,3,3,3,2))
V0<-aperm(V0)
V0[,,,,,,,,,1]<-temp[,1]
V0[,,,,,,,,,2]<-temp[,2]
V0[,,,,,,,,,3]<-temp[,3]
V0<-aperm(V0)

skip<-1
if(skip==0){
  ###WILLINGNESS TO PAY FOR POLICY REFORMS, OLD HEALTH PROBS
  newguess2<-newguess
  newguess2[c("theta_sev","eta","eta_mod","eta_sev")]<-c(-0.448,-0.185,0,0)
  param_input2<-copy(param_input)
  param_input2[c("theta_mod","theta_sev","eta","eta_single","eta_mod","eta_sev")]<-c(-0.448/2,-0.448,-0.185,-0.185,0,0)
  simdata<-Valsim_solve(as.vector(param_input2),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                        shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                        marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT)
  simdata<-cleansim(simdata,select=select,selectC=selectC)
  
  wtp_statedep<-evalwtp(paramvals=newguess2,
                       marriageprob=marriageprob, healthprob=healthprob,
                       foodstamps=foodstamps, numkids=numkids, oecd=oecd,
                       Wage=Wage, spWage=spWage, select=select,
                       selectC=selectC,
                       marriageprob_init=marriageprob_init,
                       shocks.wage = shocks.wage, shocks.spwage = shocks.spwage,
                       wagepars=wagepars, spwagepars=spwagepars,
                       onlyFT=1,evalasset=FALSE, onlySS=FALSE,onlymeanstest=FALSE,
                       single_sameallowanceprob=single_sameallowanceprob,
                       health_simpleprefs=health_simpleprefs, single_sameprefs=single_sameprefs,
                       filename="statedep")
  saveRDS(wtp_statedep,paste0("./modeloutput/wtp_statedep_",modeltype,".RDS"))
  
  rm(newguess2)
  rm(param_input2)
} else { wtp_statedep<-readRDS(paste0("./modeloutput/wtp_statedep_",modeltype,".RDS")) }

skip<-1
if(skip==0){
  ###WILLINGNESS TO PAY FOR POLICY REFORMS, FIXING MARITAL STATUS
  marriageprob_fixed<-copy(marriageprob)
  marriageprob_init_fixed<-c(0.4687035,0.5790486,0.7218947)
  for(aa in 1:3){
    for(bb in 1:3){
      marriageprob_fixed[[aa]][[bb]][,1]<-1
      marriageprob_fixed[[aa]][[bb]][,2]<-0
    }
  }

  params2<-params
  simdata<-Valsim_solve(as.vector(param_input),params2,healthprob, marriageprob_fixed,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                        shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                        marriageprob_init_fixed,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT)
  simdata<-cleansim(simdata,select=select,selectC=selectC)
  
  wtp_fixmar<-evalwtp(paramvals=newguess,
                        marriageprob=marriageprob_fixed, healthprob=healthprob,
                        foodstamps=foodstamps, numkids=numkids, oecd=oecd,
                        Wage=Wage, spWage=spWage, select=select,
                        selectC=selectC,
                        marriageprob_init=marriageprob_init_fixed,
                        shocks.wage = shocks.wage, shocks.spwage = shocks.spwage,
                        wagepars=wagepars, spwagepars=spwagepars,
                        onlyFT=1,evalasset=FALSE, onlySS=FALSE,onlymeanstest=FALSE,
                        single_sameallowanceprob=single_sameallowanceprob,
                        health_simpleprefs=health_simpleprefs, single_sameprefs=single_sameprefs,
                        filename="fixmar")
  saveRDS(wtp_fixmar,paste0("./modeloutput/wtp_fixmar_",modeltype,".RDS"))
  
  rm(newguess2)
} else { wtp_fixmar<-readRDS(paste0("./modeloutput/wtp_fixmar_",modeltype,".RDS")) }

skip<-0
if(skip==0){
  modeltype_keep<-modeltype
  numsims_keep<-copy(numsims)
  modeltype<-"twotypes"
  
  
states2<-statesetup(dir="./output",filesuffix="twotypes",
                     asset_gridsize=20,wage_gridsize=10)

  foodstamps2<-states2$foodstamps
  numkids2<-states2$numkids
#  marriageprob2<-states2$marriageprob
  marriageprob2<-marriageprob
  healthprob2<-states2$healthprob
  oecd2<-states$oecd
  
  
  #"hack" to get the two-type model to work without needing to restructur things
  #(I now have an "empty" type 1)
  marriageprob_init2<-c(0,(marriageprob_init[1]*numsims[1]+marriageprob_init[2]*numsims[2])/(numsims[1]+numsims[2]),marriageprob_init[3])
    foodstamps2[[1]]<-foodstamps2[[2]]
    numkids2[[1]]<-numkids2[[2]]
    healthprob2[[1]]<-healthprob2[[2]]
    marriageprob2[[1]]<-marriageprob2[[2]]
    oecd2[[1]]<-oecd2[[2]]

  ###WILLINGNESS TO PAY FOR POLICY REFORMS, OLD HEALTH PROBS
  numsims2<-c(0,numsims[2]+numsims[1],numsims[3])
  numsims<-copy(numsims2)
  simdata<-Valsim_solve(as.vector(param_input),params,healthprob2, marriageprob2,foodstamps2,numkids2,oecd2,aperm(Asset),aperm(Wage), aperm(spWage),
                        shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                        marriageprob_init2,numsims2,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT)
  simdata<-cleansim(simdata,select=select,selectC=selectC)
  
  wtp_notype1<-evalwtp(paramvals=newguess,
                       marriageprob=marriageprob2, healthprob=healthprob2,
                       foodstamps=foodstamps2, numkids=numkids2, oecd=oecd2,
                       Wage=Wage, spWage=spWage, select=select,
                       selectC=selectC,
                       marriageprob_init=marriageprob_init2,
                       shocks.wage = shocks.wage, shocks.spwage = shocks.spwage,
                       wagepars=wagepars, spwagepars=spwagepars,
                       onlyFT=onlyFT,evalasset=FALSE, onlySS=FALSE,onlymeanstest=FALSE,
                       single_sameallowanceprob=single_sameallowanceprob,
                       health_simpleprefs=health_simpleprefs, single_sameprefs=single_sameprefs,
                       filename="notype1")
  saveRDS(wtp_notype1,paste0("./modeloutput/wtp_notype1_",modeltype,".RDS"))
  
  modeltype<-modeltype_keep
  numsims<-copy(numsims_keep)
  rm(numsims_keep)
  rm(numsims2)
  rm(modeltype_keep)
  rm(states2,foodstampts2,numkids2,healthprob2,marriageprob2,marriageprob_init2,oecd2)

  #WTP TO BE ON DI
} else {
  modeltype_keep<-modeltype
  modeltype<-"twotypes"
  wtp_notype1<-readRDS(paste0("./modeloutput/wtp_notype1_",modeltype,".RDS")) 
  modeltype<-modeltype_keep
  rm(modeltype_keep)
}


#  wtp$lump[DI==1 & jobloss==0 & spjobloss == 0,wtp_onDI:=Assetval_solve(V0,Type,time,i,wage,spwage,asset,DI,emp,spemp,Health,spHealth,V0,Wage,spWage,Asset,policyname="getDI",wtpDI=1),by=row]
#  wtp$lump[DI==1 & jobloss==0 & spjobloss == 0,wtp_onDI_married:=Assetval_solve(V0,Type,time,i,wage,spwage,asset,DI,emp,1,Health,spHealth,V0,Wage,spWage,Asset,policyname="getDI",wtpDI=1),by=row]
#  wtp$lump[DI==1 & jobloss==0 & spjobloss == 0,wtp_onDI_single:=Assetval_solve(V0,Type,time,i,wage,spwage,asset,DI,emp,2,Health,spHealth,V0,Wage,spWage,Asset,policyname="getDI",wtpDI=1),by=row]

#plotwtp(wtp$lump,plotname="")
output<-list()

temp<-wtpchange(wtp,wtp_notype1,Type=TRUE,marital=TRUE,wtpname = "notype1",path="modeloutput/wtp_decompose")

output$summary_origcompare<-rbind(output$summary_origcompare,temp$summary_origcompare,fill=TRUE)
output$summary_typecompare<-rbind(output$summary_typecompare,temp$summary_typecompare,fill=TRUE)
output$summary_marcompare<-rbind(output$summary_marcompare,temp$summary_marcompare,fill=TRUE)
output$summary_martypecompare<-rbind(output$summary_martypecompare,temp$summary_martypecompare,fill=TRUE)

temp<-wtpchange(wtp,wtp_fixmar,Type=TRUE,marital=TRUE,wtpname = "fixmar",path="modeloutput/wtp_decompose")

output$summary_origcompare<-rbind(output$summary_origcompare,temp$summary_origcompare,fill=TRUE)
output$summary_typecompare<-rbind(output$summary_typecompare,temp$summary_typecompare,fill=TRUE)
output$summary_marcompare<-rbind(output$summary_marcompare,temp$summary_marcompare,fill=TRUE)
output$summary_martypecompare<-rbind(output$summary_martypecompare,temp$summary_martypecompare,fill=TRUE)

gc()

saveRDS(output,"./modeloutput/wtp_decompose.RDS")





  
  




